/* -*- Mode: C; tab-width: 8; indent-tabs-mode: t; c-basic-offset: 8 -*- *//**************************************************************** * * The author of this software is David M. Gay. * * Copyright (c) 1991, 2000, 2001 by Lucent Technologies. * * Permission to use, copy, modify, and distribute this software for any * purpose without fee is hereby granted, provided that this entire notice * is included in all copies of any software which is or includes a copy * or modification of this software and in all copies of the supporting * documentation for such software. * * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED * WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR LUCENT MAKES ANY * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE. * ***************************************************************//* Please send bug reports to David M. Gay (dmg at acm dot org, * with " at " changed at "@" and " dot " changed to "."). *//* On a machine with IEEE extended-precision registers, it is * necessary to specify double-precision (53-bit) rounding precision * before invoking strtod or dtoa. If the machine uses (the equivalent * of) Intel 80x87 arithmetic, the call * _control87(PC_53, MCW_PC); * does this with many compilers. Whether this or another call is * appropriate depends on the compiler; for this to work, it may be * necessary to #include "float.h" or another system-dependent header * file. *//* strtod for IEEE-, VAX-, and IBM-arithmetic machines. * * This strtod returns a nearest machine number to the input decimal * string (or sets errno to ERANGE). With IEEE arithmetic, ties are * broken by the IEEE round-even rule. Otherwise ties are broken by * biased rounding (add half and chop). * * Inspired loosely by William D. Clinger's paper "How to Read Floating * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101]. * * Modifications: * * 1. We only require IEEE, IBM, or VAX double-precision * arithmetic (not IEEE double-extended). * 2. We get by with floating-point arithmetic in a case that * Clinger missed -- when we're computing d * 10^n * for a small integer d and the integer n is not too * much larger than 22 (the maximum integer k for which * we can represent 10^k exactly), we may be able to * compute (d*10^k) * 10^(e-k) with just one roundoff. * 3. Rather than a bit-at-a-time adjustment of the binary * result in the hard case, we use floating-point * arithmetic to determine the adjustment to within * one bit; only in really hard cases do we need to * compute a second residual. * 4. Because of 3., we don't need a large table of powers of 10 * for ten-to-e (just some small tables, e.g. of 10^k * for 0 <= k <= 22). *//* * #define IEEE_8087 for IEEE-arithmetic machines where the least * significant byte has the lowest address. * #define IEEE_MC68k for IEEE-arithmetic machines where the most * significant byte has the lowest address. * #define Long int on machines with 32-bit ints and 64-bit longs. * #define IBM for IBM mainframe-style floating-point arithmetic. * #define VAX for VAX-style floating-point arithmetic (D_floating). * #define No_leftright to omit left-right logic in fast floating-point * computation of dtoa. * #define Honor_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3 * and strtod and dtoa should round accordingly. * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3 * and Honor_FLT_ROUNDS is not #defined. * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines * that use extended-precision instructions to compute rounded * products and quotients) with IBM. * #define ROUND_BIASED for IEEE-format with biased rounding. * #define Inaccurate_Divide for IEEE-format with correctly rounded * products but inaccurate quotients, e.g., for Intel i860. * #define NO_LONG_LONG on machines that do not have a "long long" * integer type (of >= 64 bits). On such machines, you can * #define Just_16 to store 16 bits per 32-bit Long when doing * high-precision integer arithmetic. Whether this speeds things * up or slows things down depends on the machine and the number * being converted. If long long is available and the name is * something other than "long long", #define Llong to be the name, * and if "unsigned Llong" does not work as an unsigned version of * Llong, #define #ULLong to be the corresponding unsigned type. * #define KR_headers for old-style C function headers. * #define Bad_float_h if your system lacks a float.h or if it does not * define some or all of DBL_DIG, DBL_MAX_10_EXP, DBL_MAX_EXP, * FLT_RADIX, FLT_ROUNDS, and DBL_MAX. * #define MALLOC your_malloc, where your_malloc(n) acts like malloc(n) * if memory is available and otherwise does something you deem * appropriate. If MALLOC is undefined, malloc will be invoked * directly -- and assumed always to succeed. Similarly, if you * want something other than the system's free() to be called to * recycle memory acquired from MALLOC, #define FREE to be the * name of the alternate routine. (Unless you #define * NO_GLOBAL_STATE and call destroydtoa, FREE or free is only * called in pathological cases, e.g., in a dtoa call after a dtoa * return in mode 3 with thousands of digits requested.) * #define Omit_Private_Memory to omit logic (added Jan. 1998) for making * memory allocations from a private pool of memory when possible. * When used, the private pool is PRIVATE_MEM bytes long: 2304 bytes, * unless #defined to be a different length. This default length * suffices to get rid of MALLOC calls except for unusual cases, * such as decimal-to-binary conversion of a very long string of * digits. The longest string dtoa can return is about 751 bytes * long. For conversions by strtod of strings of 800 digits and * all dtoa conversions in single-threaded executions with 8-byte * pointers, PRIVATE_MEM >= 7400 appears to suffice; with 4-byte * pointers, PRIVATE_MEM >= 7112 appears adequate. * #define MULTIPLE_THREADS if the system offers preemptively scheduled * multiple threads. In this case, you must provide (or suitably * #define) two locks, acquired by ACQUIRE_DTOA_LOCK(n) and freed * by FREE_DTOA_LOCK(n) for n = 0 or 1. (The second lock, accessed * in pow5mult, ensures lazy evaluation of only one copy of high * powers of 5; omitting this lock would introduce a small * probability of wasting memory, but would otherwise be harmless.) * You must also invoke freedtoa(s) to free the value s returned by * dtoa. You may do so whether or not MULTIPLE_THREADS is #defined. * #define NO_IEEE_Scale to disable new (Feb. 1997) logic in strtod that * avoids underflows on inputs whose result does not underflow. * If you #define NO_IEEE_Scale on a machine that uses IEEE-format * floating-point numbers and flushes underflows to zero rather * than implementing gradual underflow, then you must also #define * Sudden_Underflow. * #define USE_LOCALE to use the current locale's decimal_point value. * #define SET_INEXACT if IEEE arithmetic is being used and extra * computation should be done to set the inexact flag when the * result is inexact and avoid setting inexact when the result * is exact. In this case, dtoa.c must be compiled in * an environment, perhaps provided by #include "dtoa.c" in a * suitable wrapper, that defines two functions, * int get_inexact(void); * void clear_inexact(void); * such that get_inexact() returns a nonzero value if the * inexact bit is already set, and clear_inexact() sets the * inexact bit to 0. When SET_INEXACT is #defined, strtod * also does extra computations to set the underflow and overflow * flags when appropriate (i.e., when the result is tiny and * inexact or when it is a numeric value rounded to +-infinity). * #define NO_ERRNO if strtod should not assign errno = ERANGE when * the result overflows to +-Infinity or underflows to 0. * #define NO_GLOBAL_STATE to avoid defining any non-const global or * static variables. Instead the necessary state is stored in an * opaque struct, DtoaState, a pointer to which must be passed to * every entry point. Two new functions are added to the API: * DtoaState *newdtoa(void); * void destroydtoa(DtoaState *); */#ifndef Long#define Long long#endif#ifndef ULongtypedefunsignedLongULong;#endif#ifdef DEBUG#include<stdio.h>#define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);}#endif#include<stdlib.h>#include<string.h>#ifdef USE_LOCALE#include<locale.h>#endif#ifdef MALLOC#ifdef KR_headersexternchar*MALLOC();#elseexternvoid*MALLOC(size_t);#endif#else#define MALLOC malloc#endif#ifndef FREE#define FREE free#endif#ifndef Omit_Private_Memory#ifndef PRIVATE_MEM#define PRIVATE_MEM 2304#endif#define PRIVATE_mem ((PRIVATE_MEM+sizeof(double)-1)/sizeof(double))#endif#undef IEEE_Arith#undef Avoid_Underflow#ifdef IEEE_MC68k#define IEEE_Arith#endif#ifdef IEEE_8087#define IEEE_Arith#endif#include<errno.h>#ifdef Bad_float_h#ifdef IEEE_Arith#define DBL_DIG 15#define DBL_MAX_10_EXP 308#define DBL_MAX_EXP 1024#define FLT_RADIX 2#endif /*IEEE_Arith*/#ifdef IBM#define DBL_DIG 16#define DBL_MAX_10_EXP 75#define DBL_MAX_EXP 63#define FLT_RADIX 16#define DBL_MAX 7.2370055773322621e+75#endif#ifdef VAX#define DBL_DIG 16#define DBL_MAX_10_EXP 38#define DBL_MAX_EXP 127#define FLT_RADIX 2#define DBL_MAX 1.7014118346046923e+38#endif#ifndef LONG_MAX#define LONG_MAX 2147483647#endif#else /* ifndef Bad_float_h */#include<float.h>#endif /* Bad_float_h */#ifndef __MATH_H__#include<math.h>#endif#ifndef CONST#ifdef KR_headers#define CONST /* blank */#else#define CONST const#endif#endif#if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(VAX) + defined(IBM) != 1#error "Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined."#endiftypedefunion{doubled;ULongL[2];}U;#define dval(x) ((x).d)#ifdef IEEE_8087#define word0(x) ((x).L[1])#define word1(x) ((x).L[0])#else#define word0(x) ((x).L[0])#define word1(x) ((x).L[1])#endif/* The following definition of Storeinc is appropriate for MIPS processors. * An alternative that might be better on some machines is * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff) */#if defined(IEEE_8087) + defined(VAX)#define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \((unsigned short *)a)[0] = (unsigned short)c, a++)#else#define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \((unsigned short *)a)[1] = (unsigned short)c, a++)#endif/* #define P DBL_MANT_DIG *//* Ten_pmax = floor(P*log(2)/log(5)) *//* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 *//* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) *//* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */#ifdef IEEE_Arith#define Exp_shift 20#define Exp_shift1 20#define Exp_msk1 0x100000#define Exp_msk11 0x100000#define Exp_mask 0x7ff00000#define P 53#define Bias 1023#define Emin (-1022)#define Exp_1 0x3ff00000#define Exp_11 0x3ff00000#define Ebits 11#define Frac_mask 0xfffff#define Frac_mask1 0xfffff#define Ten_pmax 22#define Bletch 0x10#define Bndry_mask 0xfffff#define Bndry_mask1 0xfffff#define LSB 1#define Sign_bit 0x80000000#define Log2P 1#define Tiny0 0#define Tiny1 1#define Quick_max 14#define Int_max 14#ifndef NO_IEEE_Scale#define Avoid_Underflow#ifdef Flush_Denorm /* debugging option */#undef Sudden_Underflow#endif#endif#ifndef Flt_Rounds#ifdef FLT_ROUNDS#define Flt_Rounds FLT_ROUNDS#else#define Flt_Rounds 1#endif#endif /*Flt_Rounds*/#ifdef Honor_FLT_ROUNDS#define Rounding rounding#undef Check_FLT_ROUNDS#define Check_FLT_ROUNDS#else#define Rounding Flt_Rounds#endif#else /* ifndef IEEE_Arith */#undef Check_FLT_ROUNDS#undef Honor_FLT_ROUNDS#undef SET_INEXACT#undef Sudden_Underflow#define Sudden_Underflow#ifdef IBM#undef Flt_Rounds#define Flt_Rounds 0#define Exp_shift 24#define Exp_shift1 24#define Exp_msk1 0x1000000#define Exp_msk11 0x1000000#define Exp_mask 0x7f000000#define P 14#define Bias 65#define Exp_1 0x41000000#define Exp_11 0x41000000#define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */#define Frac_mask 0xffffff#define Frac_mask1 0xffffff#define Bletch 4#define Ten_pmax 22#define Bndry_mask 0xefffff#define Bndry_mask1 0xffffff#define LSB 1#define Sign_bit 0x80000000#define Log2P 4#define Tiny0 0x100000#define Tiny1 0#define Quick_max 14#define Int_max 15#else /* VAX */#undef Flt_Rounds#define Flt_Rounds 1#define Exp_shift 23#define Exp_shift1 7#define Exp_msk1 0x80#define Exp_msk11 0x800000#define Exp_mask 0x7f80#define P 56#define Bias 129#define Exp_1 0x40800000#define Exp_11 0x4080#define Ebits 8#define Frac_mask 0x7fffff#define Frac_mask1 0xffff007f#define Ten_pmax 24#define Bletch 2#define Bndry_mask 0xffff007f#define Bndry_mask1 0xffff007f#define LSB 0x10000#define Sign_bit 0x8000#define Log2P 1#define Tiny0 0x80#define Tiny1 0#define Quick_max 15#define Int_max 15#endif /* IBM, VAX */#endif /* IEEE_Arith */#ifndef IEEE_Arith#define ROUND_BIASED#endif#ifdef RND_PRODQUOT#define rounded_product(a,b) a = rnd_prod(a, b)#define rounded_quotient(a,b) a = rnd_quot(a, b)#ifdef KR_headersexterndoublernd_prod(),rnd_quot();#elseexterndoublernd_prod(double,double),rnd_quot(double,double);#endif#else#define rounded_product(a,b) a *= b#define rounded_quotient(a,b) a /= b#endif#define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))#define Big1 0xffffffff#ifndef Pack_32#define Pack_32#endif#ifdef KR_headers#define FFFFFFFF ((((unsigned long)0xffff)<<16)|(unsigned long)0xffff)#else#define FFFFFFFF 0xffffffffUL#endif#ifdef NO_LONG_LONG#undef ULLong#ifdef Just_16#undef Pack_32/* When Pack_32 is not defined, we store 16 bits per 32-bit Long. * This makes some inner loops simpler and sometimes saves work * during multiplications, but it often seems to make things slightly * slower. Hence the default is now to store 32 bits per Long. */#endif#else /* long long available */#ifndef Llong#define Llong long long#endif#ifndef ULLong#define ULLong unsigned Llong#endif#endif /* NO_LONG_LONG */#ifndef MULTIPLE_THREADS#define ACQUIRE_DTOA_LOCK(n) /*nothing*/#define FREE_DTOA_LOCK(n) /*nothing*/#endif#define Kmax 7structBigint{structBigint*next;intk,maxwds,sign,wds;ULongx[1];};typedefstructBigintBigint;#ifdef NO_GLOBAL_STATE#ifdef MULTIPLE_THREADS#error "cannot have both NO_GLOBAL_STATE and MULTIPLE_THREADS"#endifstructDtoaState{#define DECLARE_GLOBAL_STATE /* nothing */#else#define DECLARE_GLOBAL_STATE static#endifDECLARE_GLOBAL_STATEBigint*freelist[Kmax+1];DECLARE_GLOBAL_STATEBigint*p5s;#ifndef Omit_Private_MemoryDECLARE_GLOBAL_STATEdoubleprivate_mem[PRIVATE_mem];DECLARE_GLOBAL_STATEdouble*pmem_next#ifndef NO_GLOBAL_STATE=private_mem#endif;#endif#ifdef NO_GLOBAL_STATE};typedefstructDtoaStateDtoaState;#ifdef KR_headers#define STATE_PARAM state,#define STATE_PARAM_DECL DtoaState *state;#else#define STATE_PARAM DtoaState *state,#endif#define PASS_STATE state,#define GET_STATE(field) (state->field)staticDtoaState*newdtoa(void){DtoaState*state=(DtoaState*)MALLOC(sizeof(DtoaState));if(state){memset(state,0,sizeof(DtoaState));#ifndef Omit_Private_Memorystate->pmem_next=state->private_mem;#endif}returnstate;}staticvoiddestroydtoa#ifdef KR_headers(state)STATE_PARAM_DECL#else(DtoaState*state)#endif{inti;Bigint*v,*next;for(i=0;i<=Kmax;i++){for(v=GET_STATE(freelist)[i];v;v=next){next=v->next;#ifndef Omit_Private_Memoryif((double*)v<GET_STATE(private_mem)||(double*)v>=GET_STATE(private_mem)+PRIVATE_mem)#endifFREE((void*)v);}}#ifdef Omit_Private_MemoryBigint*p5=GET_STATE(p5s);while(p5){Bigint*tmp=p5;p5=p5->next;FREE(tmp);}#endifFREE((void*)state);}#else#define STATE_PARAM /* nothing */#define STATE_PARAM_DECL /* nothing */#define PASS_STATE /* nothing */#define GET_STATE(name) name#endifstaticBigint*Balloc#ifdef KR_headers(STATE_PARAMk)STATE_PARAM_DECLintk;#else(STATE_PARAMintk)#endif{intx;Bigint*rv;#ifndef Omit_Private_Memorysize_tlen;#endifACQUIRE_DTOA_LOCK(0);/* The k > Kmax case does not need ACQUIRE_DTOA_LOCK(0), *//* but this case seems very unlikely. */if(k<=Kmax&&(rv=GET_STATE(freelist)[k]))GET_STATE(freelist)[k]=rv->next;else{x=1<<k;#ifdef Omit_Private_Memoryrv=(Bigint*)MALLOC(sizeof(Bigint)+(x-1)*sizeof(ULong));#elselen=(sizeof(Bigint)+(x-1)*sizeof(ULong)+sizeof(double)-1)/sizeof(double);if(k<=Kmax&&GET_STATE(pmem_next)-GET_STATE(private_mem)+len<=PRIVATE_mem){rv=(Bigint*)GET_STATE(pmem_next);GET_STATE(pmem_next)+=len;}elserv=(Bigint*)MALLOC(len*sizeof(double));#endifrv->k=k;rv->maxwds=x;}FREE_DTOA_LOCK(0);rv->sign=rv->wds=0;returnrv;}staticvoidBfree#ifdef KR_headers(STATE_PARAMv)STATE_PARAM_DECLBigint*v;#else(STATE_PARAMBigint*v)#endif{if(v){if(v->k>Kmax)FREE((void*)v);else{ACQUIRE_DTOA_LOCK(0);v->next=GET_STATE(freelist)[v->k];GET_STATE(freelist)[v->k]=v;FREE_DTOA_LOCK(0);}}}#define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \y->wds*sizeof(Long) + 2*sizeof(int))staticBigint*multadd#ifdef KR_headers(STATE_PARAMb,m,a)STATE_PARAM_DECLBigint*b;intm,a;#else(STATE_PARAMBigint*b,intm,inta)/* multiply by m and add a */#endif{inti,wds;#ifdef ULLongULong*x;ULLongcarry,y;#elseULongcarry,*x,y;#ifdef Pack_32ULongxi,z;#endif#endifBigint*b1;wds=b->wds;x=b->x;i=0;carry=a;do{#ifdef ULLongy=*x*(ULLong)m+carry;carry=y>>32;*x++=(ULong)y&FFFFFFFF;#else#ifdef Pack_32xi=*x;y=(xi&0xffff)*m+carry;z=(xi>>16)*m+(y>>16);carry=z>>16;*x++=(z<<16)+(y&0xffff);#elsey=*x*m+carry;carry=y>>16;*x++=y&0xffff;#endif#endif}while(++i<wds);if(carry){if(wds>=b->maxwds){b1=Balloc(PASS_STATEb->k+1);Bcopy(b1,b);Bfree(PASS_STATEb);b=b1;}b->x[wds++]=(ULong)carry;b->wds=wds;}returnb;}staticBigint*s2b#ifdef KR_headers(STATE_PARAMs,nd0,nd,y9)STATE_PARAM_DECLCONSTchar*s;intnd0,nd;ULongy9;#else(STATE_PARAMCONSTchar*s,intnd0,intnd,ULongy9)#endif{Bigint*b;inti,k;Longx,y;x=(nd+8)/9;for(k=0,y=1;x>y;y<<=1,k++);#ifdef Pack_32b=Balloc(PASS_STATEk);b->x[0]=y9;b->wds=1;#elseb=Balloc(PASS_STATEk+1);b->x[0]=y9&0xffff;b->wds=(b->x[1]=y9>>16)?2:1;#endifi=9;if(9<nd0){s+=9;dob=multadd(PASS_STATEb,10,*s++-'0');while(++i<nd0);s++;}elses+=10;for(;i<nd;i++)b=multadd(PASS_STATEb,10,*s++-'0');returnb;}staticinthi0bits#ifdef KR_headers(x)ULongx;#else(ULongx)#endif{intk=0;if(!(x&0xffff0000)){k=16;x<<=16;}if(!(x&0xff000000)){k+=8;x<<=8;}if(!(x&0xf0000000)){k+=4;x<<=4;}if(!(x&0xc0000000)){k+=2;x<<=2;}if(!(x&0x80000000)){k++;if(!(x&0x40000000))return32;}returnk;}staticintlo0bits#ifdef KR_headers(y)ULong*y;#else(ULong*y)#endif{intk;ULongx=*y;if(x&7){if(x&1)return0;if(x&2){*y=x>>1;return1;}*y=x>>2;return2;}k=0;if(!(x&0xffff)){k=16;x>>=16;}if(!(x&0xff)){k+=8;x>>=8;}if(!(x&0xf)){k+=4;x>>=4;}if(!(x&0x3)){k+=2;x>>=2;}if(!(x&1)){k++;x>>=1;if(!x)return32;}*y=x;returnk;}staticBigint*i2b#ifdef KR_headers(STATE_PARAMi)STATE_PARAM_DECLinti;#else(STATE_PARAMinti)#endif{Bigint*b;b=Balloc(PASS_STATE1);b->x[0]=i;b->wds=1;returnb;}staticBigint*mult#ifdef KR_headers(STATE_PARAMa,b)STATE_PARAM_DECLBigint*a,*b;#else(STATE_PARAMBigint*a,Bigint*b)#endif{Bigint*c;intk,wa,wb,wc;ULong*x,*xa,*xae,*xb,*xbe,*xc,*xc0;ULongy;#ifdef ULLongULLongcarry,z;#elseULongcarry,z;#ifdef Pack_32ULongz2;#endif#endifif(a->wds<b->wds){c=a;a=b;b=c;}k=a->k;wa=a->wds;wb=b->wds;wc=wa+wb;if(wc>a->maxwds)k++;c=Balloc(PASS_STATEk);for(x=c->x,xa=x+wc;x<xa;x++)*x=0;xa=a->x;xae=xa+wa;xb=b->x;xbe=xb+wb;xc0=c->x;#ifdef ULLongfor(;xb<xbe;xc0++){if((y=*xb++)){x=xa;xc=xc0;carry=0;do{z=*x++*(ULLong)y+*xc+carry;carry=z>>32;*xc++=(ULong)z&FFFFFFFF;}while(x<xae);*xc=(ULong)carry;}}#else#ifdef Pack_32for(;xb<xbe;xb++,xc0++){if(y=*xb&0xffff){x=xa;xc=xc0;carry=0;do{z=(*x&0xffff)*y+(*xc&0xffff)+carry;carry=z>>16;z2=(*x++>>16)*y+(*xc>>16)+carry;carry=z2>>16;Storeinc(xc,z2,z);}while(x<xae);*xc=carry;}if(y=*xb>>16){x=xa;xc=xc0;carry=0;z2=*xc;do{z=(*x&0xffff)*y+(*xc>>16)+carry;carry=z>>16;Storeinc(xc,z,z2);z2=(*x++>>16)*y+(*xc&0xffff)+carry;carry=z2>>16;}while(x<xae);*xc=z2;}}#elsefor(;xb<xbe;xc0++){if(y=*xb++){x=xa;xc=xc0;carry=0;do{z=*x++*y+*xc+carry;carry=z>>16;*xc++=z&0xffff;}while(x<xae);*xc=carry;}}#endif#endiffor(xc0=c->x,xc=xc0+wc;wc>0&&!*--xc;--wc);c->wds=wc;returnc;}staticBigint*pow5mult#ifdef KR_headers(STATE_PARAMb,k)STATE_PARAM_DECLBigint*b;intk;#else(STATE_PARAMBigint*b,intk)#endif{Bigint*b1,*p5,*p51;inti;staticCONSTintp05[3]={5,25,125};if((i=k&3))b=multadd(PASS_STATEb,p05[i-1],0);if(!(k>>=2))returnb;if(!(p5=GET_STATE(p5s))){/* first time */#ifdef MULTIPLE_THREADSACQUIRE_DTOA_LOCK(1);if(!(p5=p5s)){p5=p5s=i2b(625);p5->next=0;}FREE_DTOA_LOCK(1);#elsep5=GET_STATE(p5s)=i2b(PASS_STATE625);p5->next=0;#endif}for(;;){if(k&1){b1=mult(PASS_STATEb,p5);Bfree(PASS_STATEb);b=b1;}if(!(k>>=1))break;if(!(p51=p5->next)){#ifdef MULTIPLE_THREADSACQUIRE_DTOA_LOCK(1);if(!(p51=p5->next)){p51=p5->next=mult(p5,p5);p51->next=0;}FREE_DTOA_LOCK(1);#elsep51=p5->next=mult(PASS_STATEp5,p5);p51->next=0;#endif}p5=p51;}returnb;}staticBigint*lshift#ifdef KR_headers(STATE_PARAMb,k)STATE_PARAM_DECLBigint*b;intk;#else(STATE_PARAMBigint*b,intk)#endif{inti,k1,n,n1;Bigint*b1;ULong*x,*x1,*xe,z;#ifdef Pack_32n=k>>5;#elsen=k>>4;#endifk1=b->k;n1=n+b->wds+1;for(i=b->maxwds;n1>i;i<<=1)k1++;b1=Balloc(PASS_STATEk1);x1=b1->x;for(i=0;i<n;i++)*x1++=0;x=b->x;xe=x+b->wds;#ifdef Pack_32if(k&=0x1f){k1=32-k;z=0;do{*x1++=*x<<k|z;z=*x++>>k1;}while(x<xe);if((*x1=z))++n1;}#elseif(k&=0xf){k1=16-k;z=0;do{*x1++=*x<<k&0xffff|z;z=*x++>>k1;}while(x<xe);if(*x1=z)++n1;}#endifelsedo*x1++=*x++;while(x<xe);b1->wds=n1-1;Bfree(PASS_STATEb);returnb1;}staticintcmp#ifdef KR_headers(a,b)Bigint*a,*b;#else(Bigint*a,Bigint*b)#endif{ULong*xa,*xa0,*xb,*xb0;inti,j;i=a->wds;j=b->wds;#ifdef DEBUGif(i>1&&!a->x[i-1])Bug("cmp called with a->x[a->wds-1] == 0");if(j>1&&!b->x[j-1])Bug("cmp called with b->x[b->wds-1] == 0");#endifif(i-=j)returni;xa0=a->x;xa=xa0+j;xb0=b->x;xb=xb0+j;for(;;){if(*--xa!=*--xb)return*xa<*xb?-1:1;if(xa<=xa0)break;}return0;}staticBigint*diff#ifdef KR_headers(STATE_PARAMa,b)STATE_PARAM_DECLBigint*a,*b;#else(STATE_PARAMBigint*a,Bigint*b)#endif{Bigint*c;inti,wa,wb;ULong*xa,*xae,*xb,*xbe,*xc;#ifdef ULLongULLongborrow,y;#elseULongborrow,y;#ifdef Pack_32ULongz;#endif#endifi=cmp(a,b);if(!i){c=Balloc(PASS_STATE0);c->wds=1;c->x[0]=0;returnc;}if(i<0){c=a;a=b;b=c;i=1;}elsei=0;c=Balloc(PASS_STATEa->k);c->sign=i;wa=a->wds;xa=a->x;xae=xa+wa;wb=b->wds;xb=b->x;xbe=xb+wb;xc=c->x;borrow=0;#ifdef ULLongdo{y=(ULLong)*xa++-*xb++-borrow;borrow=y>>32&(ULong)1;*xc++=(ULong)y&FFFFFFFF;}while(xb<xbe);while(xa<xae){y=*xa++-borrow;borrow=y>>32&(ULong)1;*xc++=(ULong)y&FFFFFFFF;}#else#ifdef Pack_32do{y=(*xa&0xffff)-(*xb&0xffff)-borrow;borrow=(y&0x10000)>>16;z=(*xa++>>16)-(*xb++>>16)-borrow;borrow=(z&0x10000)>>16;Storeinc(xc,z,y);}while(xb<xbe);while(xa<xae){y=(*xa&0xffff)-borrow;borrow=(y&0x10000)>>16;z=(*xa++>>16)-borrow;borrow=(z&0x10000)>>16;Storeinc(xc,z,y);}#elsedo{y=*xa++-*xb++-borrow;borrow=(y&0x10000)>>16;*xc++=y&0xffff;}while(xb<xbe);while(xa<xae){y=*xa++-borrow;borrow=(y&0x10000)>>16;*xc++=y&0xffff;}#endif#endifwhile(!*--xc)wa--;c->wds=wa;returnc;}staticdoubleulp#ifdef KR_headers(x)Ux;#else(Ux)#endif{LongL;Ua;L=(word0(x)&Exp_mask)-(P-1)*Exp_msk1;#ifndef Avoid_Underflow#ifndef Sudden_Underflowif(L>0){#endif#endif#ifdef IBML|=Exp_msk1>>4;#endifword0(a)=L;word1(a)=0;#ifndef Avoid_Underflow#ifndef Sudden_Underflow}else{L=-L>>Exp_shift;if(L<Exp_shift){word0(a)=0x80000>>L;word1(a)=0;}else{word0(a)=0;L-=Exp_shift;word1(a)=L>=31?1:1<<31-L;}}#endif#endifreturndval(a);}staticdoubleb2d#ifdef KR_headers(a,e)Bigint*a;int*e;#else(Bigint*a,int*e)#endif{ULong*xa,*xa0,w,y,z;intk;Ud;#ifdef VAXULongd0,d1;#else#define d0 word0(d)#define d1 word1(d)#endifxa0=a->x;xa=xa0+a->wds;y=*--xa;#ifdef DEBUGif(!y)Bug("zero y in b2d");#endifk=hi0bits(y);*e=32-k;#ifdef Pack_32if(k<Ebits){d0=Exp_1|y>>(Ebits-k);w=xa>xa0?*--xa:0;d1=y<<((32-Ebits)+k)|w>>(Ebits-k);gotoret_d;}z=xa>xa0?*--xa:0;if(k-=Ebits){d0=Exp_1|y<<k|z>>(32-k);y=xa>xa0?*--xa:0;d1=z<<k|y>>(32-k);}else{d0=Exp_1|y;d1=z;}#elseif(k<Ebits+16){z=xa>xa0?*--xa:0;d0=Exp_1|y<<k-Ebits|z>>Ebits+16-k;w=xa>xa0?*--xa:0;y=xa>xa0?*--xa:0;d1=z<<k+16-Ebits|w<<k-Ebits|y>>16+Ebits-k;gotoret_d;}z=xa>xa0?*--xa:0;w=xa>xa0?*--xa:0;k-=Ebits+16;d0=Exp_1|y<<k+16|z<<k|w>>16-k;y=xa>xa0?*--xa:0;d1=w<<k+16|y<<k;#endifret_d:#ifdef VAXword0(d)=d0>>16|d0<<16;word1(d)=d1>>16|d1<<16;#else#undef d0#undef d1#endifreturndval(d);}staticBigint*d2b#ifdef KR_headers(STATE_PARAMd,e,bits)STATE_PARAM_DECLUd;int*e,*bits;#else(STATE_PARAMUd,int*e,int*bits)#endif{Bigint*b;intde,k;ULong*x,y,z;#ifndef Sudden_Underflowinti;#endif#ifdef VAXULongd0,d1;d0=word0(d)>>16|word0(d)<<16;d1=word1(d)>>16|word1(d)<<16;#else#define d0 word0(d)#define d1 word1(d)#endif#ifdef Pack_32b=Balloc(PASS_STATE1);#elseb=Balloc(PASS_STATE2);#endifx=b->x;z=d0&Frac_mask;d0&=0x7fffffff;/* clear sign bit, which we ignore */#ifdef Sudden_Underflowde=(int)(d0>>Exp_shift);#ifndef IBMz|=Exp_msk11;#endif#elseif((de=(int)(d0>>Exp_shift)))z|=Exp_msk1;#endif#ifdef Pack_32if((y=d1)){if((k=lo0bits(&y))){x[0]=y|z<<(32-k);z>>=k;}elsex[0]=y;#ifndef Sudden_Underflowi=#endifb->wds=(x[1]=z)?2:1;}else{k=lo0bits(&z);x[0]=z;#ifndef Sudden_Underflowi=#endifb->wds=1;k+=32;}#elseif(y=d1){if(k=lo0bits(&y))if(k>=16){x[0]=y|z<<32-k&0xffff;x[1]=z>>k-16&0xffff;x[2]=z>>k;i=2;}else{x[0]=y&0xffff;x[1]=y>>16|z<<16-k&0xffff;x[2]=z>>k&0xffff;x[3]=z>>k+16;i=3;}else{x[0]=y&0xffff;x[1]=y>>16;x[2]=z&0xffff;x[3]=z>>16;i=3;}}else{#ifdef DEBUGif(!z)Bug("Zero passed to d2b");#endifk=lo0bits(&z);if(k>=16){x[0]=z;i=0;}else{x[0]=z&0xffff;x[1]=z>>16;i=1;}k+=32;}while(!x[i])--i;b->wds=i+1;#endif#ifndef Sudden_Underflowif(de){#endif#ifdef IBM*e=(de-Bias-(P-1)<<2)+k;*bits=4*P+8-k-hi0bits(word0(d)&Frac_mask);#else*e=de-Bias-(P-1)+k;*bits=P-k;#endif#ifndef Sudden_Underflow}else{*e=de-Bias-(P-1)+1+k;#ifdef Pack_32*bits=32*i-hi0bits(x[i-1]);#else*bits=(i+2)*16-hi0bits(x[i]);#endif}#endifreturnb;}#undef d0#undef d1staticdoubleratio#ifdef KR_headers(a,b)Bigint*a,*b;#else(Bigint*a,Bigint*b)#endif{Uda,db;intk,ka,kb;dval(da)=b2d(a,&ka);dval(db)=b2d(b,&kb);#ifdef Pack_32k=ka-kb+32*(a->wds-b->wds);#elsek=ka-kb+16*(a->wds-b->wds);#endif#ifdef IBMif(k>0){word0(da)+=(k>>2)*Exp_msk1;if(k&=3)dval(da)*=1<<k;}else{k=-k;word0(db)+=(k>>2)*Exp_msk1;if(k&=3)dval(db)*=1<<k;}#elseif(k>0)word0(da)+=k*Exp_msk1;else{k=-k;word0(db)+=k*Exp_msk1;}#endifreturndval(da)/dval(db);}staticCONSTdoubletens[]={1e0,1e1,1e2,1e3,1e4,1e5,1e6,1e7,1e8,1e9,1e10,1e11,1e12,1e13,1e14,1e15,1e16,1e17,1e18,1e19,1e20,1e21,1e22#ifdef VAX,1e23,1e24#endif};staticCONSTdouble#ifdef IEEE_Arithbigtens[]={1e16,1e32,1e64,1e128,1e256};staticCONSTdoubletinytens[]={1e-16,1e-32,1e-64,1e-128,#ifdef Avoid_Underflow9007199254740992.*9007199254740992.e-256/* = 2^106 * 1e-53 */#else1e-256#endif};/* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow *//* flag unnecessarily. It leads to a song and dance at the end of strtod. */#define Scale_Bit 0x10#define n_bigtens 5#else#ifdef IBMbigtens[]={1e16,1e32,1e64};staticCONSTdoubletinytens[]={1e-16,1e-32,1e-64};#define n_bigtens 3#elsebigtens[]={1e16,1e32};staticCONSTdoubletinytens[]={1e-16,1e-32};#define n_bigtens 2#endif#endifstaticdouble_strtod#ifdef KR_headers(STATE_PARAMs00,se)STATE_PARAM_DECLCONSTchar*s00;char**se;#else(STATE_PARAMCONSTchar*s00,char**se)#endif{#ifdef Avoid_Underflowintscale;#endifintbb2,bb5,bbe,bd2,bd5,bbbits,bs2,c,dsign,e,e1,esign,i,j,k,nd,nd0,nf,nz,nz0,sign;CONSTchar*s,*s0,*s1;doubleaadj,adj;Uaadj1,rv,rv0;LongL;ULongy,z;Bigint*bb,*bb1,*bd,*bd0,*bs,*delta;#ifdef SET_INEXACTintinexact,oldinexact;#endif#ifdef Honor_FLT_ROUNDSintrounding;#endif#ifdef USE_LOCALECONSTchar*s2;#endif#ifdef __GNUC__delta=bb=bd=bs=0;#endifsign=nz0=nz=0;dval(rv)=0.;for(s=s00;;s++)switch(*s){case'-':sign=1;/* no break */case'+':if(*++s)gotobreak2;/* no break */case0:gotoret0;case'\t':case'\n':case'\v':case'\f':case'\r':case' ':continue;default:gotobreak2;}break2:if(*s=='0'){nz0=1;while(*++s=='0');if(!*s)gotoret;}s0=s;y=z=0;for(nd=nf=0;(c=*s)>='0'&&c<='9';nd++,s++)if(nd<9)y=10*y+c-'0';elseif(nd<16)z=10*z+c-'0';nd0=nd;#ifdef USE_LOCALEs1=localeconv()->decimal_point;if(c==*s1){c='.';if(*++s1){s2=s;for(;;){if(*++s2!=*s1){c=0;break;}if(!*++s1){s=s2;break;}}}}#endifif(c=='.'){c=*++s;if(!nd){for(;c=='0';c=*++s)nz++;if(c>'0'&&c<='9'){s0=s;nf+=nz;nz=0;gotohave_dig;}gotodig_done;}for(;c>='0'&&c<='9';c=*++s){have_dig:nz++;if(c-='0'){nf+=nz;for(i=1;i<nz;i++)if(nd++<9)y*=10;elseif(nd<=DBL_DIG+1)z*=10;if(nd++<9)y=10*y+c;elseif(nd<=DBL_DIG+1)z=10*z+c;nz=0;}}}dig_done:e=0;if(c=='e'||c=='E'){if(!nd&&!nz&&!nz0){gotoret0;}s00=s;esign=0;switch(c=*++s){case'-':esign=1;case'+':c=*++s;}if(c>='0'&&c<='9'){while(c=='0')c=*++s;if(c>'0'&&c<='9'){L=c-'0';s1=s;while((c=*++s)>='0'&&c<='9')L=10*L+c-'0';if(s-s1>8||L>19999)/* Avoid confusion from exponents * so large that e might overflow. */e=19999;/* safe for 16 bit ints */elsee=(int)L;if(esign)e=-e;}elsee=0;}elses=s00;}if(!nd){if(!nz&&!nz0){ret0:s=s00;sign=0;}gotoret;}e1=e-=nf;/* Now we have nd0 digits, starting at s0, followed by a * decimal point, followed by nd-nd0 digits. The number we're * after is the integer represented by those digits times * 10**e */if(!nd0)nd0=nd;k=nd<DBL_DIG+1?nd:DBL_DIG+1;dval(rv)=y;if(k>9){#ifdef SET_INEXACTif(k>DBL_DIG)oldinexact=get_inexact();#endifdval(rv)=tens[k-9]*dval(rv)+z;}bd0=0;if(nd<=DBL_DIG#ifndef RND_PRODQUOT#ifndef Honor_FLT_ROUNDS&&Flt_Rounds==1#endif#endif){if(!e)gotoret;if(e>0){if(e<=Ten_pmax){#ifdef VAXgotovax_ovfl_check;#else#ifdef Honor_FLT_ROUNDS/* round correctly FLT_ROUNDS = 2 or 3 */if(sign){rv=-rv;sign=0;}#endif/* rv = */rounded_product(dval(rv),tens[e]);gotoret;#endif}i=DBL_DIG-nd;if(e<=Ten_pmax+i){/* A fancier test would sometimes let us do * this for larger i values. */#ifdef Honor_FLT_ROUNDS/* round correctly FLT_ROUNDS = 2 or 3 */if(sign){rv=-rv;sign=0;}#endife-=i;dval(rv)*=tens[i];#ifdef VAX/* VAX exponent range is so narrow we must * worry about overflow here... */vax_ovfl_check:word0(rv)-=P*Exp_msk1;/* rv = */rounded_product(dval(rv),tens[e]);if((word0(rv)&Exp_mask)>Exp_msk1*(DBL_MAX_EXP+Bias-1-P))gotoovfl;word0(rv)+=P*Exp_msk1;#else/* rv = */rounded_product(dval(rv),tens[e]);#endifgotoret;}}#ifndef Inaccurate_Divideelseif(e>=-Ten_pmax){#ifdef Honor_FLT_ROUNDS/* round correctly FLT_ROUNDS = 2 or 3 */if(sign){rv=-rv;sign=0;}#endif/* rv = */rounded_quotient(dval(rv),tens[-e]);gotoret;}#endif}e1+=nd-k;#ifdef IEEE_Arith#ifdef SET_INEXACTinexact=1;if(k<=DBL_DIG)oldinexact=get_inexact();#endif#ifdef Avoid_Underflowscale=0;#endif#ifdef Honor_FLT_ROUNDSif((rounding=Flt_Rounds)>=2){if(sign)rounding=rounding==2?0:2;elseif(rounding!=2)rounding=0;}#endif#endif /*IEEE_Arith*//* Get starting approximation = rv * 10**e1 */if(e1>0){if((i=e1&15))dval(rv)*=tens[i];if(e1&=~15){if(e1>DBL_MAX_10_EXP){ovfl:#ifndef NO_ERRNOerrno=ERANGE;#endif/* Can't trust HUGE_VAL */#ifdef IEEE_Arith#ifdef Honor_FLT_ROUNDSswitch(rounding){case0:/* toward 0 */case3:/* toward -infinity */word0(rv)=Big0;word1(rv)=Big1;break;default:word0(rv)=Exp_mask;word1(rv)=0;}#else /*Honor_FLT_ROUNDS*/word0(rv)=Exp_mask;word1(rv)=0;#endif /*Honor_FLT_ROUNDS*/#ifdef SET_INEXACT/* set overflow bit */dval(rv0)=1e300;dval(rv0)*=dval(rv0);#endif#else /*IEEE_Arith*/word0(rv)=Big0;word1(rv)=Big1;#endif /*IEEE_Arith*/if(bd0)gotoretfree;gotoret;}e1>>=4;for(j=0;e1>1;j++,e1>>=1)if(e1&1)dval(rv)*=bigtens[j];/* The last multiplication could overflow. */word0(rv)-=P*Exp_msk1;dval(rv)*=bigtens[j];if((z=word0(rv)&Exp_mask)>Exp_msk1*(DBL_MAX_EXP+Bias-P))gotoovfl;if(z>Exp_msk1*(DBL_MAX_EXP+Bias-1-P)){/* set to largest number *//* (Can't trust DBL_MAX) */word0(rv)=Big0;word1(rv)=Big1;}elseword0(rv)+=P*Exp_msk1;}}elseif(e1<0){e1=-e1;if((i=e1&15))dval(rv)/=tens[i];if(e1>>=4){if(e1>=1<<n_bigtens)gotoundfl;#ifdef Avoid_Underflowif(e1&Scale_Bit)scale=2*P;for(j=0;e1>0;j++,e1>>=1)if(e1&1)dval(rv)*=tinytens[j];if(scale&&(j=2*P+1-((word0(rv)&Exp_mask)>>Exp_shift))>0){/* scaled rv is denormal; zap j low bits */if(j>=32){word1(rv)=0;if(j>=53)word0(rv)=(P+2)*Exp_msk1;elseword0(rv)&=0xffffffff<<(j-32);}elseword1(rv)&=0xffffffff<<j;}#elsefor(j=0;e1>1;j++,e1>>=1)if(e1&1)dval(rv)*=tinytens[j];/* The last multiplication could underflow. */dval(rv0)=dval(rv);dval(rv)*=tinytens[j];if(!dval(rv)){dval(rv)=2.*dval(rv0);dval(rv)*=tinytens[j];#endifif(!dval(rv)){undfl:dval(rv)=0.;#ifndef NO_ERRNOerrno=ERANGE;#endifif(bd0)gotoretfree;gotoret;}#ifndef Avoid_Underflowword0(rv)=Tiny0;word1(rv)=Tiny1;/* The refinement below will clean * this approximation up. */}#endif}}/* Now the hard part -- adjusting rv to the correct value.*//* Put digits into bd: true value = bd * 10^e */bd0=s2b(PASS_STATEs0,nd0,nd,y);for(;;){bd=Balloc(PASS_STATEbd0->k);Bcopy(bd,bd0);bb=d2b(PASS_STATErv,&bbe,&bbbits);/* rv = bb * 2^bbe */bs=i2b(PASS_STATE1);if(e>=0){bb2=bb5=0;bd2=bd5=e;}else{bb2=bb5=-e;bd2=bd5=0;}if(bbe>=0)bb2+=bbe;elsebd2-=bbe;bs2=bb2;#ifdef Honor_FLT_ROUNDSif(rounding!=1)bs2++;#endif#ifdef Avoid_Underflowj=bbe-scale;i=j+bbbits-1;/* logb(rv) */if(i<Emin)/* denormal */j+=P-Emin;elsej=P+1-bbbits;#else /*Avoid_Underflow*/#ifdef Sudden_Underflow#ifdef IBMj=1+4*P-3-bbbits+((bbe+bbbits-1)&3);#elsej=P+1-bbbits;#endif#else /*Sudden_Underflow*/j=bbe;i=j+bbbits-1;/* logb(rv) */if(i<Emin)/* denormal */j+=P-Emin;elsej=P+1-bbbits;#endif /*Sudden_Underflow*/#endif /*Avoid_Underflow*/bb2+=j;bd2+=j;#ifdef Avoid_Underflowbd2+=scale;#endifi=bb2<bd2?bb2:bd2;if(i>bs2)i=bs2;if(i>0){bb2-=i;bd2-=i;bs2-=i;}if(bb5>0){bs=pow5mult(PASS_STATEbs,bb5);bb1=mult(PASS_STATEbs,bb);Bfree(PASS_STATEbb);bb=bb1;}if(bb2>0)bb=lshift(PASS_STATEbb,bb2);if(bd5>0)bd=pow5mult(PASS_STATEbd,bd5);if(bd2>0)bd=lshift(PASS_STATEbd,bd2);if(bs2>0)bs=lshift(PASS_STATEbs,bs2);delta=diff(PASS_STATEbb,bd);dsign=delta->sign;delta->sign=0;i=cmp(delta,bs);#ifdef Honor_FLT_ROUNDSif(rounding!=1){if(i<0){/* Error is less than an ulp */if(!delta->x[0]&&delta->wds<=1){/* exact */#ifdef SET_INEXACTinexact=0;#endifbreak;}if(rounding){if(dsign){adj=1.;gotoapply_adj;}}elseif(!dsign){adj=-1.;if(!word1(rv)&&!(word0(rv)&Frac_mask)){y=word0(rv)&Exp_mask;#ifdef Avoid_Underflowif(!scale||y>2*P*Exp_msk1)#elseif(y)#endif{delta=lshift(PASS_STATEdelta,Log2P);if(cmp(delta,bs)<=0)adj=-0.5;}}apply_adj:#ifdef Avoid_Underflowif(scale&&(y=word0(rv)&Exp_mask)<=2*P*Exp_msk1)word0(adj)+=(2*P+1)*Exp_msk1-y;#else#ifdef Sudden_Underflowif((word0(rv)&Exp_mask)<=P*Exp_msk1){word0(rv)+=P*Exp_msk1;dval(rv)+=adj*ulp(rv);word0(rv)-=P*Exp_msk1;}else#endif /*Sudden_Underflow*/#endif /*Avoid_Underflow*/dval(rv)+=adj*ulp(rv);}break;}adj=ratio(delta,bs);if(adj<1.)adj=1.;if(adj<=0x7ffffffe){/* adj = rounding ? ceil(adj) : floor(adj); */y=adj;if(y!=adj){if(!((rounding>>1)^dsign))y++;adj=y;}}#ifdef Avoid_Underflowif(scale&&(y=word0(rv)&Exp_mask)<=2*P*Exp_msk1)word0(adj)+=(2*P+1)*Exp_msk1-y;#else#ifdef Sudden_Underflowif((word0(rv)&Exp_mask)<=P*Exp_msk1){word0(rv)+=P*Exp_msk1;adj*=ulp(rv);if(dsign)dval(rv)+=adj;elsedval(rv)-=adj;word0(rv)-=P*Exp_msk1;gotocont;}#endif /*Sudden_Underflow*/#endif /*Avoid_Underflow*/adj*=ulp(rv);if(dsign)dval(rv)+=adj;elsedval(rv)-=adj;gotocont;}#endif /*Honor_FLT_ROUNDS*/if(i<0){/* Error is less than half an ulp -- check for * special case of mantissa a power of two. */if(dsign||word1(rv)||word0(rv)&Bndry_mask#ifdef IEEE_Arith#ifdef Avoid_Underflow||(word0(rv)&Exp_mask)<=(2*P+1)*Exp_msk1#else||(word0(rv)&Exp_mask)<=Exp_msk1#endif#endif){#ifdef SET_INEXACTif(!delta->x[0]&&delta->wds<=1)inexact=0;#endifbreak;}if(!delta->x[0]&&delta->wds<=1){/* exact result */#ifdef SET_INEXACTinexact=0;#endifbreak;}delta=lshift(PASS_STATEdelta,Log2P);if(cmp(delta,bs)>0)gotodrop_down;break;}if(i==0){/* exactly half-way between */if(dsign){if((word0(rv)&Bndry_mask1)==Bndry_mask1&&word1(rv)==(#ifdef Avoid_Underflow(scale&&(y=word0(rv)&Exp_mask)<=2*P*Exp_msk1)?(0xffffffff&(0xffffffff<<(2*P+1-(y>>Exp_shift)))):#endif0xffffffff)){/*boundary case -- increment exponent*/word0(rv)=(word0(rv)&Exp_mask)+Exp_msk1#ifdef IBM|Exp_msk1>>4#endif;word1(rv)=0;#ifdef Avoid_Underflowdsign=0;#endifbreak;}}elseif(!(word0(rv)&Bndry_mask)&&!word1(rv)){drop_down:/* boundary case -- decrement exponent */#ifdef Sudden_Underflow /*{{*/L=word0(rv)&Exp_mask;#ifdef IBMif(L<Exp_msk1)#else#ifdef Avoid_Underflowif(L<=(scale?(2*P+1)*Exp_msk1:Exp_msk1))#elseif(L<=Exp_msk1)#endif /*Avoid_Underflow*/#endif /*IBM*/gotoundfl;L-=Exp_msk1;#else /*Sudden_Underflow}{*/#ifdef Avoid_Underflowif(scale){L=word0(rv)&Exp_mask;if(L<=(2*P+1)*Exp_msk1){if(L>(P+2)*Exp_msk1)/* round even ==> *//* accept rv */break;/* rv = smallest denormal */gotoundfl;}}#endif /*Avoid_Underflow*/L=(word0(rv)&Exp_mask)-Exp_msk1;#endif /*Sudden_Underflow}}*/word0(rv)=L|Bndry_mask1;word1(rv)=0xffffffff;#ifdef IBMgotocont;#elsebreak;#endif}#ifndef ROUND_BIASEDif(!(word1(rv)&LSB))break;#endifif(dsign)dval(rv)+=ulp(rv);#ifndef ROUND_BIASEDelse{dval(rv)-=ulp(rv);#ifndef Sudden_Underflowif(!dval(rv))gotoundfl;#endif}#ifdef Avoid_Underflowdsign=1-dsign;#endif#endifbreak;}if((aadj=ratio(delta,bs))<=2.){if(dsign)aadj=dval(aadj1)=1.;elseif(word1(rv)||word0(rv)&Bndry_mask){#ifndef Sudden_Underflowif(word1(rv)==Tiny1&&!word0(rv))gotoundfl;#endifaadj=1.;dval(aadj1)=-1.;}else{/* special case -- power of FLT_RADIX to be *//* rounded down... */if(aadj<2./FLT_RADIX)aadj=1./FLT_RADIX;elseaadj*=0.5;dval(aadj1)=-aadj;}}else{aadj*=0.5;dval(aadj1)=dsign?aadj:-aadj;#ifdef Check_FLT_ROUNDSswitch(Rounding){case2:/* towards +infinity */dval(aadj1)-=0.5;break;case0:/* towards 0 */case3:/* towards -infinity */dval(aadj1)+=0.5;}#elseif(Flt_Rounds==0)dval(aadj1)+=0.5;#endif /*Check_FLT_ROUNDS*/}y=word0(rv)&Exp_mask;/* Check for overflow */if(y==Exp_msk1*(DBL_MAX_EXP+Bias-1)){dval(rv0)=dval(rv);word0(rv)-=P*Exp_msk1;adj=dval(aadj1)*ulp(rv);dval(rv)+=adj;if((word0(rv)&Exp_mask)>=Exp_msk1*(DBL_MAX_EXP+Bias-P)){if(word0(rv0)==Big0&&word1(rv0)==Big1)gotoovfl;word0(rv)=Big0;word1(rv)=Big1;gotocont;}elseword0(rv)+=P*Exp_msk1;}else{#ifdef Avoid_Underflowif(scale&&y<=2*P*Exp_msk1){if(aadj<=0x7fffffff){if((z=(ULong)aadj)<=0)z=1;aadj=z;dval(aadj1)=dsign?aadj:-aadj;}word0(aadj1)+=(2*P+1)*Exp_msk1-y;}adj=dval(aadj1)*ulp(rv);dval(rv)+=adj;#else#ifdef Sudden_Underflowif((word0(rv)&Exp_mask)<=P*Exp_msk1){dval(rv0)=dval(rv);word0(rv)+=P*Exp_msk1;adj=dval(aadj1)*ulp(rv);dval(rv)+=adj;#ifdef IBMif((word0(rv)&Exp_mask)<P*Exp_msk1)#elseif((word0(rv)&Exp_mask)<=P*Exp_msk1)#endif{if(word0(rv0)==Tiny0&&word1(rv0)==Tiny1)gotoundfl;word0(rv)=Tiny0;word1(rv)=Tiny1;gotocont;}elseword0(rv)-=P*Exp_msk1;}else{adj=dval(aadj1)*ulp(rv);dval(rv)+=adj;}#else /*Sudden_Underflow*//* Compute adj so that the IEEE rounding rules will * correctly round rv + adj in some half-way cases. * If rv * ulp(rv) is denormalized (i.e., * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid * trouble from bits lost to denormalization; * example: 1.2e-307 . */if(y<=(P-1)*Exp_msk1&&aadj>1.){dval(aadj1)=(double)(int)(aadj+0.5);if(!dsign)dval(aadj1)=-dval(aadj1);}adj=dval(aadj1)*ulp(rv);dval(rv)+=adj;#endif /*Sudden_Underflow*/#endif /*Avoid_Underflow*/}z=word0(rv)&Exp_mask;#ifndef SET_INEXACT#ifdef Avoid_Underflowif(!scale)#endifif(y==z){/* Can we stop now? */L=(Long)aadj;aadj-=L;/* The tolerances below are conservative. */if(dsign||word1(rv)||word0(rv)&Bndry_mask){if(aadj<.4999999||aadj>.5000001)break;}elseif(aadj<.4999999/FLT_RADIX)break;}#endifcont:Bfree(PASS_STATEbb);Bfree(PASS_STATEbd);Bfree(PASS_STATEbs);Bfree(PASS_STATEdelta);}#ifdef SET_INEXACTif(inexact){if(!oldinexact){word0(rv0)=Exp_1+(70<<Exp_shift);word1(rv0)=0;dval(rv0)+=1.;}}elseif(!oldinexact)clear_inexact();#endif#ifdef Avoid_Underflowif(scale){word0(rv0)=Exp_1-2*P*Exp_msk1;word1(rv0)=0;dval(rv)*=dval(rv0);#ifndef NO_ERRNO/* try to avoid the bug of testing an 8087 register value */if(word0(rv)==0&&word1(rv)==0)errno=ERANGE;#endif}#endif /* Avoid_Underflow */#ifdef SET_INEXACTif(inexact&&!(word0(rv)&Exp_mask)){/* set underflow bit */dval(rv0)=1e-300;dval(rv0)*=dval(rv0);}#endifretfree:Bfree(PASS_STATEbb);Bfree(PASS_STATEbd);Bfree(PASS_STATEbs);Bfree(PASS_STATEbd0);Bfree(PASS_STATEdelta);ret:if(se)*se=(char*)s;returnsign?-dval(rv):dval(rv);}staticintquorem#ifdef KR_headers(b,S)Bigint*b,*S;#else(Bigint*b,Bigint*S)#endif{intn;ULong*bx,*bxe,q,*sx,*sxe;#ifdef ULLongULLongborrow,carry,y,ys;#elseULongborrow,carry,y,ys;#ifdef Pack_32ULongsi,z,zs;#endif#endifn=S->wds;#ifdef DEBUG/*debug*/if(b->wds>n)/*debug*/Bug("oversize b in quorem");#endifif(b->wds<n)return0;sx=S->x;sxe=sx+--n;bx=b->x;bxe=bx+n;q=*bxe/(*sxe+1);/* ensure q <= true quotient */#ifdef DEBUG/*debug*/if(q>9)/*debug*/Bug("oversized quotient in quorem");#endifif(q){borrow=0;carry=0;do{#ifdef ULLongys=*sx++*(ULLong)q+carry;carry=ys>>32;y=*bx-(ys&FFFFFFFF)-borrow;borrow=y>>32&(ULong)1;*bx++=(ULong)y&FFFFFFFF;#else#ifdef Pack_32si=*sx++;ys=(si&0xffff)*q+carry;zs=(si>>16)*q+(ys>>16);carry=zs>>16;y=(*bx&0xffff)-(ys&0xffff)-borrow;borrow=(y&0x10000)>>16;z=(*bx>>16)-(zs&0xffff)-borrow;borrow=(z&0x10000)>>16;Storeinc(bx,z,y);#elseys=*sx++*q+carry;carry=ys>>16;y=*bx-(ys&0xffff)-borrow;borrow=(y&0x10000)>>16;*bx++=y&0xffff;#endif#endif}while(sx<=sxe);if(!*bxe){bx=b->x;while(--bxe>bx&&!*bxe)--n;b->wds=n;}}if(cmp(b,S)>=0){q++;borrow=0;carry=0;bx=b->x;sx=S->x;do{#ifdef ULLongys=*sx+++carry;carry=ys>>32;y=*bx-(ys&FFFFFFFF)-borrow;borrow=y>>32&(ULong)1;*bx++=(ULong)y&FFFFFFFF;#else#ifdef Pack_32si=*sx++;ys=(si&0xffff)+carry;zs=(si>>16)+(ys>>16);carry=zs>>16;y=(*bx&0xffff)-(ys&0xffff)-borrow;borrow=(y&0x10000)>>16;z=(*bx>>16)-(zs&0xffff)-borrow;borrow=(z&0x10000)>>16;Storeinc(bx,z,y);#elseys=*sx+++carry;carry=ys>>16;y=*bx-(ys&0xffff)-borrow;borrow=(y&0x10000)>>16;*bx++=y&0xffff;#endif#endif}while(sx<=sxe);bx=b->x;bxe=bx+n;if(!*bxe){while(--bxe>bx&&!*bxe)--n;b->wds=n;}}returnq;}#if !defined(MULTIPLE_THREADS) && !defined(NO_GLOBAL_STATE)#define USE_DTOA_RESULT 1staticchar*dtoa_result;#endifstaticchar*#ifdef KR_headersrv_alloc(STATE_PARAMi)STATE_PARAM_DECLinti;#elserv_alloc(STATE_PARAMinti)#endif{intj,k,*r;j=sizeof(ULong);for(k=0;sizeof(Bigint)-sizeof(ULong)-sizeof(int)+j<=(unsigned)i;j<<=1)k++;r=(int*)Balloc(PASS_STATEk);*r=k;return#ifdef USE_DTOA_RESULTdtoa_result=#endif(char*)(r+1);}staticchar*#ifdef KR_headersnrv_alloc(STATE_PARAMs,rve,n)STATE_PARAM_DECLchar*s,**rve;intn;#elsenrv_alloc(STATE_PARAMCONSTchar*s,char**rve,intn)#endif{char*rv,*t;t=rv=rv_alloc(PASS_STATEn);while((*t=*s++))t++;if(rve)*rve=t;returnrv;}/* freedtoa(s) must be used to free values s returned by dtoa * when MULTIPLE_THREADS is #defined. It should be used in all cases, * but for consistency with earlier versions of dtoa, it is optional * when MULTIPLE_THREADS is not defined. */staticvoid#ifdef KR_headersfreedtoa(STATE_PARAMs)STATE_PARAM_DECLchar*s;#elsefreedtoa(STATE_PARAMchar*s)#endif{Bigint*b=(Bigint*)((int*)s-1);b->maxwds=1<<(b->k=*(int*)b);Bfree(PASS_STATEb);#ifdef USE_DTOA_RESULTif(s==dtoa_result)dtoa_result=0;#endif}/* dtoa for IEEE arithmetic (dmg): convert double to ASCII string. * * Inspired by "How to Print Floating-Point Numbers Accurately" by * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126]. * * Modifications: * 1. Rather than iterating, we use a simple numeric overestimate * to determine k = floor(log10(d)). We scale relevant * quantities using O(log2(k)) rather than O(k) multiplications. * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't * try to generate digits strictly left to right. Instead, we * compute with fewer bits and propagate the carry if necessary * when rounding the final digit up. This is often faster. * 3. Under the assumption that input will be rounded nearest, * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22. * That is, we allow equality in stopping tests when the * round-nearest rule will give the same floating-point value * as would satisfaction of the stopping test with strict * inequality. * 4. We remove common factors of powers of 2 from relevant * quantities. * 5. When converting floating-point integers less than 1e16, * we use floating-point arithmetic rather than resorting * to multiple-precision integers. * 6. When asked to produce fewer than 15 digits, we first try * to get by with floating-point arithmetic; we resort to * multiple-precision integer arithmetic only if we cannot * guarantee that the floating-point calculation has given * the correctly rounded result. For k requested digits and * "uniformly" distributed input, the probability is * something like 10^(k-15) that we must resort to the Long * calculation. */staticchar*dtoa#ifdef KR_headers(STATE_PARAMd,mode,ndigits,decpt,sign,rve)STATE_PARAM_DECLUd;intmode,ndigits,*decpt,*sign;char**rve;#else(STATE_PARAMUd,intmode,intndigits,int*decpt,int*sign,char**rve)#endif{/* Arguments ndigits, decpt, sign are similar to those of ecvt and fcvt; trailing zeros are suppressed from the returned string. If not null, *rve is set to point to the end of the return value. If d is +-Infinity or NaN, then *decpt is set to 9999. mode: 0 ==> shortest string that yields d when read in and rounded to nearest. 1 ==> like 0, but with Steele & White stopping rule; e.g. with IEEE P754 arithmetic , mode 0 gives 1e23 whereas mode 1 gives 9.999999999999999e22. 2 ==> max(1,ndigits) significant digits. This gives a return value similar to that of ecvt, except that trailing zeros are suppressed. 3 ==> through ndigits past the decimal point. This gives a return value similar to that from fcvt, except that trailing zeros are suppressed, and ndigits can be negative. 4,5 ==> similar to 2 and 3, respectively, but (in round-nearest mode) with the tests of mode 0 to possibly return a shorter string that rounds to d. With IEEE arithmetic and compilation with -DHonor_FLT_ROUNDS, modes 4 and 5 behave the same as modes 2 and 3 when FLT_ROUNDS != 1. 6-9 ==> Debugging modes similar to mode - 4: don't try fast floating-point estimate (if applicable). Values of mode other than 0-9 are treated as mode 0. Sufficient space is allocated to the return value to hold the suppressed trailing zeros. */intbbits,b2,b5,be,dig,i,ieps,ilim,ilim0,ilim1,j,j1,k,k0,k_check,leftright,m2,m5,s2,s5,spec_case,try_quick;LongL;#ifndef Sudden_Underflowintdenorm;ULongx;#endifBigint*b,*b1,*delta,*mlo,*mhi,*S;Ud2,eps;doubleds;char*s,*s0;#ifdef Honor_FLT_ROUNDSintrounding;#endif#ifdef SET_INEXACTintinexact,oldinexact;#endif#ifdef __GNUC__ilim=ilim1=0;mlo=NULL;#endif#ifdef USE_DTOA_RESULTif(dtoa_result){freedtoa(PASS_STATEdtoa_result);dtoa_result=0;}#endifif(word0(d)&Sign_bit){/* set sign for everything, including 0's and NaNs */*sign=1;word0(d)&=~Sign_bit;/* clear sign bit */}else*sign=0;#if defined(IEEE_Arith) + defined(VAX)#ifdef IEEE_Arithif((word0(d)&Exp_mask)==Exp_mask)#elseif(word0(d)==0x8000)#endif{/* Infinity or NaN */*decpt=9999;#ifdef IEEE_Arithif(!word1(d)&&!(word0(d)&0xfffff))returnnrv_alloc(PASS_STATE"Infinity",rve,8);#endifreturnnrv_alloc(PASS_STATE"NaN",rve,3);}#endif#ifdef IBMdval(d)+=0;/* normalize */#endifif(!dval(d)){*decpt=1;returnnrv_alloc(PASS_STATE"0",rve,1);}#ifdef SET_INEXACTtry_quick=oldinexact=get_inexact();inexact=1;#endif#ifdef Honor_FLT_ROUNDSif((rounding=Flt_Rounds)>=2){if(*sign)rounding=rounding==2?0:2;elseif(rounding!=2)rounding=0;}#endifb=d2b(PASS_STATEd,&be,&bbits);#ifdef Sudden_Underflowi=(int)(word0(d)>>Exp_shift1&(Exp_mask>>Exp_shift1));#elseif((i=(int)(word0(d)>>Exp_shift1&(Exp_mask>>Exp_shift1)))){#endifdval(d2)=dval(d);word0(d2)&=Frac_mask1;word0(d2)|=Exp_11;#ifdef IBMif(j=11-hi0bits(word0(d2)&Frac_mask))dval(d2)/=1<<j;#endif/* log(x) ~=~ log(1.5) + (x-1.5)/1.5 * log10(x) = log(x) / log(10) * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10)) * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2) * * This suggests computing an approximation k to log10(d) by * * k = (i - Bias)*0.301029995663981 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 ); * * We want k to be too large rather than too small. * The error in the first-order Taylor series approximation * is in our favor, so we just round up the constant enough * to compensate for any error in the multiplication of * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077, * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14, * adding 1e-13 to the constant term more than suffices. * Hence we adjust the constant term to 0.1760912590558. * (We could get a more accurate k by invoking log10, * but this is probably not worthwhile.) */i-=Bias;#ifdef IBMi<<=2;i+=j;#endif#ifndef Sudden_Underflowdenorm=0;}else{/* d is denormalized */i=bbits+be+(Bias+(P-1)-1);x=i>32?word0(d)<<(64-i)|word1(d)>>(i-32):word1(d)<<(32-i);dval(d2)=x;word0(d2)-=31*Exp_msk1;/* adjust exponent */i-=(Bias+(P-1)-1)+1;denorm=1;}#endifds=(dval(d2)-1.5)*0.289529654602168+0.1760912590558+i*0.301029995663981;k=(int)ds;if(ds<0.&&ds!=k)k--;/* want k = floor(ds) */k_check=1;if(k>=0&&k<=Ten_pmax){if(dval(d)<tens[k])k--;k_check=0;}j=bbits-i-1;if(j>=0){b2=0;s2=j;}else{b2=-j;s2=0;}if(k>=0){b5=0;s5=k;s2+=k;}else{b2-=k;b5=-k;s5=0;}if(mode<0||mode>9)mode=0;#ifndef SET_INEXACT#ifdef Check_FLT_ROUNDStry_quick=Rounding==1;#elsetry_quick=1;#endif#endif /*SET_INEXACT*/if(mode>5){mode-=4;try_quick=0;}leftright=1;switch(mode){case0:case1:ilim=ilim1=-1;i=18;ndigits=0;break;case2:leftright=0;/* no break */case4:if(ndigits<=0)ndigits=1;ilim=ilim1=i=ndigits;break;case3:leftright=0;/* no break */case5:i=ndigits+k+1;ilim=i;ilim1=i-1;if(i<=0)i=1;}s=s0=rv_alloc(PASS_STATEi);#ifdef Honor_FLT_ROUNDSif(mode>1&&rounding!=1)leftright=0;#endifif(ilim>=0&&ilim<=Quick_max&&try_quick){/* Try to get by with floating-point arithmetic. */i=0;dval(d2)=dval(d);k0=k;ilim0=ilim;ieps=2;/* conservative */if(k>0){ds=tens[k&0xf];j=k>>4;if(j&Bletch){/* prevent overflows */j&=Bletch-1;dval(d)/=bigtens[n_bigtens-1];ieps++;}for(;j;j>>=1,i++)if(j&1){ieps++;ds*=bigtens[i];}dval(d)/=ds;}elseif((j1=-k)){dval(d)*=tens[j1&0xf];for(j=j1>>4;j;j>>=1,i++)if(j&1){ieps++;dval(d)*=bigtens[i];}}if(k_check&&dval(d)<1.&&ilim>0){if(ilim1<=0)gotofast_failed;ilim=ilim1;k--;dval(d)*=10.;ieps++;}dval(eps)=ieps*dval(d)+7.;word0(eps)-=(P-1)*Exp_msk1;if(ilim==0){S=mhi=0;dval(d)-=5.;if(dval(d)>dval(eps))gotoone_digit;if(dval(d)<-dval(eps))gotono_digits;gotofast_failed;}#ifndef No_leftrightif(leftright){/* Use Steele & White method of only * generating digits needed. */dval(eps)=0.5/tens[ilim-1]-dval(eps);for(i=0;;){L=(ULong)dval(d);dval(d)-=L;*s++='0'+(int)L;if(dval(d)<dval(eps))gotoret1;if(1.-dval(d)<dval(eps))gotobump_up;if(++i>=ilim)break;dval(eps)*=10.;dval(d)*=10.;}}else{#endif/* Generate ilim digits, then fix them up. */dval(eps)*=tens[ilim-1];for(i=1;;i++,dval(d)*=10.){L=(Long)(dval(d));if(!(dval(d)-=L))ilim=i;*s++='0'+(int)L;if(i==ilim){if(dval(d)>0.5+dval(eps))gotobump_up;elseif(dval(d)<0.5-dval(eps)){while(*--s=='0');s++;gotoret1;}break;}}#ifndef No_leftright}#endiffast_failed:s=s0;dval(d)=dval(d2);k=k0;ilim=ilim0;}/* Do we have a "small" integer? */if(be>=0&&k<=Int_max){/* Yes. */ds=tens[k];if(ndigits<0&&ilim<=0){S=mhi=0;if(ilim<0||dval(d)<5*ds)gotono_digits;gotoone_digit;}for(i=1;;i++,dval(d)*=10.){L=(Long)(dval(d)/ds);dval(d)-=L*ds;#ifdef Check_FLT_ROUNDS/* If FLT_ROUNDS == 2, L will usually be high by 1 */if(dval(d)<0){L--;dval(d)+=ds;}#endif*s++='0'+(int)L;if(!dval(d)){#ifdef SET_INEXACTinexact=0;#endifbreak;}if(i==ilim){#ifdef Honor_FLT_ROUNDSif(mode>1)switch(rounding){case0:gotoret1;case2:gotobump_up;}#endifdval(d)+=dval(d);if(dval(d)>ds||(dval(d)==ds&&L&1)){bump_up:while(*--s=='9')if(s==s0){k++;*s='0';break;}++*s++;}break;}}gotoret1;}m2=b2;m5=b5;mhi=mlo=0;if(leftright){i=#ifndef Sudden_Underflowdenorm?be+(Bias+(P-1)-1+1):#endif#ifdef IBM1+4*P-3-bbits+((bbits+be-1)&3);#else1+P-bbits;#endifb2+=i;s2+=i;mhi=i2b(PASS_STATE1);}if(m2>0&&s2>0){i=m2<s2?m2:s2;b2-=i;m2-=i;s2-=i;}if(b5>0){if(leftright){if(m5>0){mhi=pow5mult(PASS_STATEmhi,m5);b1=mult(PASS_STATEmhi,b);Bfree(PASS_STATEb);b=b1;}if((j=b5-m5))b=pow5mult(PASS_STATEb,j);}elseb=pow5mult(PASS_STATEb,b5);}S=i2b(PASS_STATE1);if(s5>0)S=pow5mult(PASS_STATES,s5);/* Check for special case that d is a normalized power of 2. */spec_case=0;if((mode<2||leftright)#ifdef Honor_FLT_ROUNDS&&rounding==1#endif){if(!word1(d)&&!(word0(d)&Bndry_mask)#ifndef Sudden_Underflow&&word0(d)&(Exp_mask&~Exp_msk1)#endif){/* The special case */b2+=Log2P;s2+=Log2P;spec_case=1;}}/* Arrange for convenient computation of quotients: * shift left if necessary so divisor has 4 leading 0 bits. * * Perhaps we should just compute leading 28 bits of S once * and for all and pass them and a shift to quorem, so it * can do shifts and ors to compute the numerator for q. */#ifdef Pack_32if((i=((s5?32-hi0bits(S->x[S->wds-1]):1)+s2)&0x1f))i=32-i;#elseif(i=((s5?32-hi0bits(S->x[S->wds-1]):1)+s2)&0xf)i=16-i;#endifif(i>4){i-=4;b2+=i;m2+=i;s2+=i;}elseif(i<4){i+=28;b2+=i;m2+=i;s2+=i;}if(b2>0)b=lshift(PASS_STATEb,b2);if(s2>0)S=lshift(PASS_STATES,s2);if(k_check){if(cmp(b,S)<0){k--;b=multadd(PASS_STATEb,10,0);/* we botched the k estimate */if(leftright)mhi=multadd(PASS_STATEmhi,10,0);ilim=ilim1;}}if(ilim<=0&&(mode==3||mode==5)){if(ilim<0||cmp(b,S=multadd(PASS_STATES,5,0))<0){/* no digits, fcvt style */no_digits:/* MOZILLA CHANGE: Always return a non-empty string. */*s++='0';k=0;gotoret;}one_digit:*s++='1';k++;gotoret;}if(leftright){if(m2>0)mhi=lshift(PASS_STATEmhi,m2);/* Compute mlo -- check for special case * that d is a normalized power of 2. */mlo=mhi;if(spec_case){mhi=Balloc(PASS_STATEmhi->k);Bcopy(mhi,mlo);mhi=lshift(PASS_STATEmhi,Log2P);}for(i=1;;i++){dig=quorem(b,S)+'0';/* Do we yet have the shortest decimal string * that will round to d? */j=cmp(b,mlo);delta=diff(PASS_STATES,mhi);j1=delta->sign?1:cmp(b,delta);Bfree(PASS_STATEdelta);#ifndef ROUND_BIASEDif(j1==0&&mode!=1&&!(word1(d)&1)#ifdef Honor_FLT_ROUNDS&&rounding>=1#endif){if(dig=='9')gotoround_9_up;if(j>0)dig++;#ifdef SET_INEXACTelseif(!b->x[0]&&b->wds<=1)inexact=0;#endif*s++=dig;gotoret;}#endifif(j<0||(j==0&&mode!=1#ifndef ROUND_BIASED&&!(word1(d)&1)#endif)){if(!b->x[0]&&b->wds<=1){#ifdef SET_INEXACTinexact=0;#endifgotoaccept_dig;}#ifdef Honor_FLT_ROUNDSif(mode>1)switch(rounding){case0:gotoaccept_dig;case2:gotokeep_dig;}#endif /*Honor_FLT_ROUNDS*/if(j1>0){b=lshift(PASS_STATEb,1);j1=cmp(b,S);if((j1>0||(j1==0&&dig&1))&&dig++=='9')gotoround_9_up;}accept_dig:*s++=dig;gotoret;}if(j1>0){#ifdef Honor_FLT_ROUNDSif(!rounding)gotoaccept_dig;#endifif(dig=='9'){/* possible if i == 1 */round_9_up:*s++='9';gotoroundoff;}*s++=dig+1;gotoret;}#ifdef Honor_FLT_ROUNDSkeep_dig:#endif*s++=dig;if(i==ilim)break;b=multadd(PASS_STATEb,10,0);if(mlo==mhi)mlo=mhi=multadd(PASS_STATEmhi,10,0);else{mlo=multadd(PASS_STATEmlo,10,0);mhi=multadd(PASS_STATEmhi,10,0);}}}elsefor(i=1;;i++){*s++=dig=quorem(b,S)+'0';if(!b->x[0]&&b->wds<=1){#ifdef SET_INEXACTinexact=0;#endifgotoret;}if(i>=ilim)break;b=multadd(PASS_STATEb,10,0);}/* Round off last digit */#ifdef Honor_FLT_ROUNDSswitch(rounding){case0:gototrimzeros;case2:gotoroundoff;}#endifb=lshift(PASS_STATEb,1);j=cmp(b,S);if(j>=0){/* ECMA compatible rounding needed by Spidermonkey */roundoff:while(*--s=='9')if(s==s0){k++;*s++='1';gotoret;}++*s++;}else{#ifdef Honor_FLT_ROUNDStrimzeros:#endifwhile(*--s=='0');s++;}ret:Bfree(PASS_STATES);if(mhi){if(mlo&&mlo!=mhi)Bfree(PASS_STATEmlo);Bfree(PASS_STATEmhi);}ret1:#ifdef SET_INEXACTif(inexact){if(!oldinexact){word0(d)=Exp_1+(70<<Exp_shift);word1(d)=0;dval(d)+=1.;}}elseif(!oldinexact)clear_inexact();#endifBfree(PASS_STATEb);*s=0;*decpt=k+1;if(rve)*rve=s;returns0;}#undef CONST